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With the forecast increase in air traffic demand over the next decades, it is imper- 
ative to develop tools to provide traffic flow managers with the information required 
to support decision making. In particular, decision-support tools for traffic flow man- 
agement should aid in limiting controller workload and complexity, while supporting 
increases in air traffic throughput. Indeed, the growth of air transportation is con- 
ditioned on the ability to maintain acceptable safety. While mainy decision-support 
tools exist for short-term traffic planning, few have addressed the strategic needs for 
medium- and long-term planning for time horizons greater than 30 minutes. This pa- 
per seeks to address this gap through the introduction of 3D aircraft proximity maps 
that evaluate the future probability of presence of at least one or two aircraft at any 
given point of the airspace. Three types of proximity maps are presented: (i) pres- 
ence maps that indicate the local density of traffic, (ii) conflict maps that determine 
locations of potential conflicts and their corresponding probabilities, and (iii) outlier 
proximity maps that show the probability of conflict due to aircraft not belonging to 
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dominant traffic patterns. These maps provide traffic flow managers with information 
relating to the complexity and difficulty of managing an airspace. The intended pur- 
pose of the maps is to anticipate how aircraft flows interact, and how outliers impact 
the dominant traffic flow for a given time period. This formulation is able to predict 
which "critical" regions may be subject to conflicts between aircraft, thereby requiring 
careful monitoring and additional effort to manage the airspace. These probabilities 
are computed using a generative aircraft flow model. Time-varying flow character- 
istics, such as geometrical configuration, speed, and probability density function of 
aircraft spatial distribution within the flow, are determined from archived Enhanced 
Traffic Management System data, using a tailored clustering algorithm. Aircraft not 
belonging to flows are identifled as outliers. 

I. Introduction 

ACCORDING to the predicted air traffic growth, current airspace capacity Hmits. i.e. the 
maximum number of aircraft allowed in an airspace, will be increasingly exceeded throughout the 
National Airspace System (NAS) in the coming decades Currently, the capacity of an airspace is 
measured according to the peak number of aircraft allowed in the considered sector within any one 
minute period 0. When the predetermined thresholds are exceeded, congestion alerts are triggered. 
In practice, capacity thresholds are set by the ability of controllers to handle traffic and to maintain 
manageable workloads for worst-case planning over a wide range of traffic patterns [3 . Hence, 
current airspace capacity thresholds are independent of any changes in the traffic configurations 
that may occur due to weather or other extraneous traffic perturbations. Such a ill-defined method 
for establishing and measuring airspace capacity often leads to under-utilization of the airspace 
under nominal conditions. Or conversely, when convective weather is present, over-utilization is not 
uncommon. To accommodate high levels of air traffic throughput and to effectively manage airspace 
utilization while maintaining safety, intrinsic descriptors of traffic and the airspace are required to 
support airspace management in the redesign or adjustment of flow configurations. These measures 
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are often referred to as "complexity measures", since they aim at providing an image of the difficulty 
or complexity of managing the considered airspace and traffic configuration. 

There has been a significant volume of research related to the estimation of air traffic control 
complexity, and its application to determine airspace capacity. One method for defining and es- 
tablishing air traffic control complexity is dynamic density [4, 5J. This metric is a weighted sum 
of the local density of aircraft; the number of heading, altitude, and speed changes; and a variety 
of other factors. The weighting factors are determined and validated through human-in-the-loop 
experiments. Dynamic density makes use of instantaneous aircraft configurations, and therefore 
requires accurate prediction of traffic to hold any significant forecasting value. The impact of un- 
certainties on aircraft trajectories on predicting dynamic density has been studied in [5j[6], and was 
demonstrated to cause short-comings in forecasting complexity. Another intrinsic measure of air 
traffic control complexity, proposed in [7_, analyzes the stability of a nonlinear dynamical system 
corresponding to the considered traffic pattern. This analysis method is attractive for its ability 
to seamlessly account for aircraft position uncertainties. An alternative approach toward managing 
complexity is presented in [8j. The authors propose a dynamic sectorization of the airspace based 
on air traffic controller workload, which leads to the generation of another kind of complexity map. 

Recently, significant effort has been placed into understanding the probability of conflict and 
the associated resolution-commands required to ensure separation. For example, several tools are 
presented in [3WTT] to determine the probability of conflict for a given geometric traffic configuration, 
with and without aircraft position uncertainty. In |12| . an "input-output" approach is proposed that 
generates complexity maps based on the conflict-resolution commands required to accommodate 
disturbances such as the entrance of an exogenous aircraft into the airspace. In [13', the authors use 
aircraft position and intent information to compute the probability of presence of the aircraft at any 
point in the airspace within a given time frame. In |T¥lll5| . the authors propose a complexity metric 
that is a function of the local average velocity vector, the density of aircraft and the proximity of 
an aircraft to the considered point of the airspace, in order to evaluate the solutions computed by 
their traffic flow management (TFM) algorithm. 

Nevertheless, a key piece of information forming the foundation of the previously cited works is 
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knowledge of the initial aircraft positions and trajectories ("known" with or without uncertainty). 
Such information, while useful for short-term air traffic management, e.g. planning a safe trajectory 
for an intruder aircraft through the airspace (e.g. [131) j is not well-suited for mid-term or long-term 
air traffic management, that is, for time horizons greater than 30 minutes. Indeed, from a Traffic 
Flow Management perspective, strategic planning requires knowledge of air traffic patterns and 
flows characteristics (e.g. geometrical conflguration, flight plans, distances between two consecutive 
aircraft), and not the current position and intent information of individual aircraft. As explained 
in considering aircraft flows enables the computation of more predictable and perturbation- 
resistant estimates than considering individual aircraft. Representing air traffic as flows provides 
a more meaningful and higher-level view of air traffic to traffic flow managers. In fact, such a 
representation is also consistent with an air traffic controllers understanding of a sector. As explained 
in [17L I18| , one such way air traffic controllers abstract sectors is according to dominant traffic flow 
patterns and traffic flow interactions. Few works have developed TFM tools with this viewpoint in 
mind. The authors of fl6 provide a list of flow features to describe flow pattern, such as the number 
of flows, the major flow and its size, the number of crossing flows, etc. In [,19._, these criteria allow 
the authors to estimate and predict sector capacity as a function of the traffic flow pattern, and 
to study the impact of severe weather. While the list of flow features established in |TB] provides 
a notion of the traffic complexity (given as a scalar number), more information may be needed 
by traffic flow managers to identify regions of high complexity. In [20 , a simplified approach to 
generate complexity maps, taking into account the influence of automated an conflict-resolution 
algorithm and the flows characteristics, is proposed. However, only two flows are considered, and 
the extension of the results of |20| to a realistic airspace (i.e. multiple intersections with multiple 
flows) is not obvious and may be computationally expensive. 

Building upon the preliminary work in [21j , this paper proposes a new kind of air traffic com- 
plexity maps aimed at supporting traffic flow management. These maps, referred to as "(aircraft) 
proximity maps", provide an estimation of probability of proximity of one or more aircraft to a 
region of airspace. For simplicity, the expression probability of proximity is used throughout the 
paper instead of the correct (and more rigorous) expression probability an aircraft is within a given 
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neighborhood of a point P in the airspace. Three kinds of proximity maps are developed, each with 
a different probability metric. They are as follows: 

• Presence maps depict the probability of presence of at least one aircraft from a flow in a given 
volume of the airspace. These maps enable the identification of the main flows of aircraft and 
the regions of high aircraft density. 

• Conflict maps provide the probability of presence of at least two aircraft from two different 
flows in a given volume of the airspace. These maps highlight regions of airspace with high 
probability of conflict or close proximity between aircraft, thereby indicating regions of the 
airspace that require additional attention from the air traffic controllers. 

• Outlier proximity maps give the probability of presence of an aircraft from the dominant traffic 



and an outlier aircrafl (deflned in Section HE) in a considered volume of the airspace. These 



maps aid traffic flow managers in evaluating the impact of aircraft outliers on the dominant 
air traffic flows. 

Regions of the maps with high values for the probability of presence, conflict, or outlier proximity 
are interpreted to be "critical" regions that require monitoring. 

The proposed proximity maps take into account flow characteristics in a manner consistent 
with TFM (i.e. aggregate flow descriptions such as average miles-in-trail) , thereby making them an 
intuitive and effective strategic decision-support tool. To enable the creation of the presence and 
proximity maps, a consistence model of airspace is required. However, because a describing airspace 
is quite complex, we introduce a dimensionally-reduced model. The dimensionally-reduced model 
serves as a generative model that is consistent with the tools we are proposing. That is, the initial 
modeling provides a framework in which a parametric model can be established. For this purpose, 
the clustering algorithm introduced in [22_ is used to process archived Enhanced Traffic Management 
System (ETMS) data to generate a realistic air traffic flow model, i.e. the generative model. The 
resulting model accounts for the aircraft spatial distributions and aircraft arrival distributions for 
each flow. Ultimately, the methodology for generating the traffic model provides a data-driven 
framework on which proximity maps are based. The process, starting from archived data to the 
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Fig. 1 Work-process flow from historical data to the generation of proximity maps. 



generation of proximity maps, is illustrated in Fig. [T] with an example of clustered trajectories and 
presence maps. 

The usage of proximity maps extends beyond displaying critical information. They may be 
used for decision making in traffic flow management. Examples include dynamically adjusting 
flow characteristics; redefining or creating new flight routes; or visualizing the influence of rogue 
aircraft (similar to [T^) or severe weather (as in P^) located in a specific airspace region. The 
tools presented in this article may also be utilized as decision support tools in conjunction with the 
Flow-Based Route Planning method used in |14l [T5] . or other such tools. 

The remainder of the paper is organized as follows. Section |IT] presents a methodology for 
developing a traffic model through the extraction of air traffic flows and outliers from archived 
flight data. The resulting model serves as a parametric generative model for computation of the 
proximity maps. In the next section, an empirical analysis of the clustering is completed. The 
analysis is completed to support assumptions and simplifying calculations used in the proximity 
maps. Section [TV] introduces the notion of the "residual distance" of an aircraft and describes how 
presence, conflict and outlier proximity maps are generated from the flow and outlier models. In 
Section [Vj example proximity maps are generated for a representative region of airspace. Finally, 
in Section [VI[ our conclusions are presented. 



II. Airspace Modeling: Transformation of flight data to flow and outlier models 

This section introduces a methodology to build mathematical models of aircraft operations 
within an airspace derived from ETMS data. The resulting model serves as a generative parametric 
model to create the proposed proximity maps. To precede, useful definitions are first introduced. 
Next, the ETMS data set is described, followed by a brief discussion of the clustering algorithm 
used to generate the traffic models. The clustering algorithm is then applied to an typical airspace, 
and analysis of the resulting traffic flow and outlier model is provided. 

A. Deflnitions 

Consider a 3D airspace ^ S M'^, whose arbitrarily shaped boundaries are fixed and independent 
of the altitude. The coordinates of a point P — (x, y, z) of the airspace are expressed in an East- 
North- Up frame TZq, where the x,y axes are expressed in nautical miles (NM) and the altitude z 
is expressed in feet (ft); a simple latitude-longitude to flat earth projection is used, with the origin 
of the airspace is located at the center of the airspace. For the scope and scale of the airspace 
(approximately 400 NM in the largest dimension), the error introduced by this approximation is 
insignificant in relation to the topological descriptors need to understand generate the proximity 
maps. As customary. Flight Levels (FL) are used for the z-coordinate. Flight levels are discretized 
at 100 ft increment, such that FL 320 corresponds to an altitude of 32,000 ft. 

The following terms, used throughout the article, some of which are depicted in Fig. [2] are 
defined as follows: 

• A trajectory is a sequence of I points defined by a set of 3D coordinates. All trajectories are 
resampled to have same number of points, I ~ 8. Trajectories are spatially sampled, instead 
of temporally sampled, in order to allow for identification of spatially similar trajectories. 

• A cluster, Ci, is a labeling of rrii trajectories that are identified as "similar" by the clustering 
algorithm. 

• A track % is the centroid of the cluster Ci, i.e the trajectory "center of mass" of all the 
trajectories in the cluster. The track 7i can be interpreted as the representative trajectory for 
cluster Ci. 
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• The window of cluster Ci is located at the k^^ point along the associated track %. The 
window Wi is vertically orientated with its face aligned according to the heading of Ti at the 
fc"^ point towards the (fc + 1)*^ point. The dimensions of each window are determined by the 
extrema of the fc*^ points corresponding to each trajectory within the cluster labeling. 

• At each window W^, the local coordinate frame (a|, b^, c|) located at point fc**^ point on track 
T/'. The vertical axis, c|, is aligned with the global z coordinate axis, while a| align with the 
track Ti, such that it points towards W^~^^. 

• A 60a; Bf,k = 1 ... Z — 1, is the polyhedron defined by the convex hull of two consecutive 
windows, W,* and W^+\ 

• A tube is the succession of all the boxes of a given cluster. 

• For every window , the sampled lateral distribution of aircraft, f^^k, is computed from the 
orthogonal projections of the k^^ point of the trajectories in Ci onto the axis of window 

• For every window Wf , the sampled vertical distribution of aircraft, f^k is computed from the 
altitudes of the k^^ point of the trajectories in Cj along the cf axis. 

• The aircraft arrival rate parameter X{t), is the expected number of aircraft entering the 
airspace region during a given time period. The arrival rate for each cluster Ci, A,(t), is 
defined similarly. For the remainder of the paper, the considered time period is 15 minutes 
long. 

• A flow J^i is a 3-dimensional tube that contains the trajectories of cluster C^. The tube is 
defined by a succession of I windows , = 1 . . .1. Associated with a flow are its aircraft 
speed distribution, aircraft arrival rate, and its aircraft lateral and vertical distributions. 

B. ETMS Data Set 

The data used to construct the airspace traffic model is taken from the Enhanced Traflic Manage- 
ment System (ETMS) data set. Cleveland center is selected for the study because of its significance 
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Fig. 2 A 3D representation of an ascending flow with its corresponding tube. 

to the NAS. Indeed, a large share of the traffic to and from the North-East of the United States flies 
through the center. Furthermore, as noted in [23 , Cleveland sector is "one of the most congested 
and delay-prone Centers in the continental United States," and centrally to located to a number of 
high traffic regions. 

The ETMS data includes aircraft radar tracks (longitude, latitude, altitude) taken at approx- 
imately 1 minute intervals. During the 123 days (May to August 2005) covered by the data, all 
526,840 aircraft trajectories with at least one point over FL 250 are considered; this subset of the 
data covers the majority of en route aircraft. The dates 08/04/2005, 08/24/2005, and 08/26/2005 
are not included in the clustering analysis, as the corresponding data was not available. ETMS data 
is notorious for its inconsistencies in altitudes, therefore a simple filtering algorithm on reported 
altitudes was used. After this initial filtering process, aircraft with inconsistencies in reported posi- 
tions and altitudes are excluded. Flight data is deemed to be inconsistent if it exhibits physically 
unrealizable trajectories, e.g. rapid changes in altitude or position. After filtering out inconsistent 
trajectories, a "clean" data-set of 338,060 trajectories remains. Statistics not requiring the altitude 
or the position information, such as the aircraft arrival rate into the airspace are computed on the 
entire en route data set, while the clustering of trajectories is run on the "clean" subset. Figure [3] 
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presents a randomly sampled density plot of 110,000 "clean" trajectories over Cleveland center (only 
110,000 flights are plotted because of memory limitations). 

Major air traffic flows are clearly identifiable as the darkest regions of the airspace; many 
corresponding to defined jetways. The appearance of the gray shades indicates that some aircraft 
fly outside of the major flows. The clustering algorithm presented in the next subsection is used to 
identify the major flows, thereby establishing to dominant traffic pattern. 

Density plot for 110,000 trajectories 




.200' ' ' ' ' ' ' ' ' 
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X-axis, NM 

Fig. 3 Density plot of a subset of the "clean" trajectories. 

C. Trajectory Clustering 

The 338,060 trajectories are clustered using a slightly modified version of the algorithm presented 
in The algorithm is defined according to the following steps: 

1. Clean and format the trajectories. Inconsistent trajectories are excluded. The remaining 
trajectories are spatially resampled according to a linear interpolation with a fixed number of 
equally space points. The linear approximation is sufficient as most aircraft fly along straight 
trajectories in this center. 

2. Augment the dimensionality of the trajectory data by including additional features such as 
heading, polar coordinates, etc. 
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3. Apply hierarchical clustering. Organize and separate the trajectories by altitude and attitude 
to create separate data sets. Level trajectories are split by flight level; climbing trajectories 
are sorted by destination (or airspace region exit) flight level; and descending trajectories are 
split by initial flight level. 

4. Normalize each feature and concatenate the data into a single row vector for each flight. 
Each hierarchical grouping - deflned by altitude and attitude - has a corresponding matrix of 
data, including interpolated trajectory points and augmented data. Each row of the matrices 
corresponds to a flight, and each column corresponds to a feature. 

5. Apply principal component analysis ^24j on the previously described matrices and keep the 
projections onto the first 5 principal components, thereby reducing the dimensionality of the 



6. Cluster the values of the projections using a density-based clustering algorithm (DB- 



7. Obtain clusters of trajectories and outliers for each altitude and attitude category. Outliers 
are trajectories that are not assigned to a cluster. 

Here, a key extension to the work presented in is to split the dataset into subsets, by 
altitude and attitude (climb, descent, level flight), to greater fldelity in modeling. This initial 
categorization results in 136 subsets. As noted, the trajectories for each category are clustered 
using the presented algorithm. The number of clusters and intra-cluster similarity can be adjusted 
by tuning input parameters to DBSCAN (the minimum number of points required to start a cluster, 
and the minimum distance between those points). With such a large sample set of trajectories, a 
two-pass approach is taken towards the clustering. The flrst round identifles both large and small 
clusters. A second application of DBSCAN is applied to large clusters with more discriminating 
parameters. 

Following categorization and clustering, approximatively 80% of the trajectories are grouped 
into 685 clusters, or flows. Figure |4] presents a 2D and 3D view of the centroids of the clusters. 
Blue lines represent westbound level traffic, yellow lines represent eastbound level traffic, green lines 
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Fig. 4 Centroids for all traffic flow clusters. 

represent descending traffic, and red lines represent ascending traffic. Eastbound traffic (yellow lines) 
appear to be solely passing through the center of the center. The reason is that Eastbound traffic is 
widespread over the entire center and no clear major flow is identiflable. Therefore, for each flight 
level, most of the Eastbound trajectories belong to the same cluster which is very wide. In Figure 
|4(b)[ major airports of Cleveland center - Cleveland Hopkins airport (CLE), Detroit Metropolitan 
Wayne County Airport (DTW) and Pittsburgh International Airport (PIT) - are clearly identiflable 
by the clusters corresponding to descending air traffic. Given the greater variability in takeoff 
trajectories, in comparison to landing trajectories, few ascending clusters from airports are identifled. 
Primary reasons for variability include the diverse range in aircraft performance and pilot decisions 
characteristics. 



D. Clustering Analysis 

Applying clustering to the flight data identifles numerous flows over a broad range of trajectory 
types. To better understand the structure of air traffic in relation to the clustering algorithm, results 
are analyzed over a broad range of parameters to determine when flow models of the airspace may 
lose validity. Of particular interest is the proportion of outlier aircraft, that is, aircraft trajectories 
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Fig. 5 Distribution of clustered nad outlier flights by flight level. 

that are not assigned a cluster labeling. 

Figure [5] presents the distribution of clustered and outlier trajectories by flight level. "Leveled at 
FL" indicates that the altitude is constant for the corresponding aircraft trajectory; "Ascending to 
FL" indicates the destination flight level for the aircraft trajectory within the center; and "Descending 
from FL" indicates the initial flight level of the aircraft before starting its descent. For the considered 
flight levels over 80% of trajectories belong to a cluster. This result supports the hypothesis that 
flows exist, and can be meaningful towards describing air traffic. 

Further analysis is completed to determine the fraction of trajectories clustered according to 
time of day. and date. Figure [6] illustrates the distribution of clustered and outlier trajectories by 
time of the day. According to the flgure. the fraction of outliers is relatively constant. Meanwhile, 
the number of trajectories clustered varies throughout the day. with the maximum number of clusters 
being established during peak hours, i.e. between 11 a.m. and 2 a.m. GMT (5 a.m. and 8 p.m. 
local time). Figure |7] presents the number of trajectories identified as part of a cluster, or identified 
as an outlier for each of the 123 days (1 = May 1, 2005). 

The results presented here demonstrate that the clustering remains consistent over a broad 
range of parameterizations (time of day, date, altitude/attitude). Therefore, clustering yields a 
model that can be utilized for subsequent complexity analysis. 
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Fig. 7 Analysis of clustered and outliers flights by day of traffic. 



E. Outliers 

While the vast majority of the traffic (approximatively 80%) are clustered into flows, 20% are 
designated to be outliers. Outliers are trajectories that are deemed "too different" from trajectories 
assigned flow labels by the clustering algorithm. "Too different" is a relative measure that is reflected 
by parameters of the clustering algorithm. In order to approach a complete model of the airspace, it 
is vital that non-clustered trajectories are accounted for. Outlier aircraft are particularly important 
to assessing air traffic control complexity because they act as perturbations to an otherwise static 
system. As such, complete traffic models require inclusion of outlier aircraft that are not clustered 
according to standard traffic flows. The outlier model is given in the next section, with the complete 
parameterized air traffic model. 
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III. Creating a Parameterized Generative Model 

The previous section provided a methodology for clustering aircraft trajectories into flows and for 
determining which aircraft can be considered as outlier (non-clustered aircraft). Using the clusters 
as a basis, the model is further analyzed to create a parameterized generative model. Conceptually, 
the parameterized generative model provides an input to generate proximity complexity maps. 
Parameterizing the airspace traffic model in this fashion allows for tools to test 'what-if ' scenarios 
or to create proximity maps without re-clustering and analyzing the air traffic data. To parameterize 
traffic, the aircraft arrival process and the speed distribution for each flow are analyzed and modeled, 
as well as the outliers trajectories. 



A. Spatial Flow Distribution 

The trajectories labeled to a particular flow by the clustering algorithm do not follow a singular 
trajectory, hence the introduction of trajectory windows and flow tubes. In fact, some flows have a 
large width, up to several hundred miles for Eastbound flight levels with a large throughput. 



Figure 8(a) presents an example cluster, with corresponding trajectories, windows, and tube. 
The vertical and lateral aircraft distributions for a speciflc window (resp. for each window W*^) are 
illustrated lateral Figure [8(c) The vertical and horizontal distributions appear to be decorrelated. 
For the generative model, the decorrelation of vertical and lateral distribution are assumed. To 
support this hypothesis. Figure |9] shows a histogram of correlation values between the vertical 
and lateral distributions, for each of the 8 windows for the 685 flows. The distribution of the 
correlation values between the vertical and lateral distributions is narrow, and centered in 0. For 
90% of all cluster windows, the absolute value of the correlation is less than 0.31. This allows the 
parameterized generative model to use the approximation that the vertical and lateral distribution 
of aircraft trajectories are independent. 



B. Speed distribution 

Within a given flow, aircraft fly at different speeds. Factors for such variation include differences 
in aircraft and engine type, aircraft weight, weather, wind, and pilot preferences. Analysis of the 
flows was conducted to determine a distribution to best model aircraft speeds. Ultimately, the t- 
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(b) Distribution analysis at a window 
Fig. 8 Lateral and vertical distributions of aircraft positions at a window within a flow. 

location scale probability distribution was selected. The probability distribution, , for the aircraft 
speeds in flow J^i is 



(1) 



where ^ is the location parameter, cr > is the scale parameter, > is the shape parameter, and 
r is the standard gamma probability distribution. If x has a t-scale distribution with parameters 



fj,,a and v, then has a Student's i-distribution with v degrees of freedom. Figure 10 presents 



an example of an empirical probability density function for aircraft speeds within a major flow, as 
well as the t-location scale fit for the empirical distribution. 
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Fig. 9 Correlation between lateral and vertical distributions of aircraft positions at the win- 
dows. 
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Fig. 10 Distribution of aircraft speeds for the busiest flow. 



C. Aircraft arrival process 

The clustering method provides a spatial representation of the flows. To obtain a complete traflic 
model it is necessary that the trafliic rate iirto each flow is appropriately described from historical 
data. Due to the sparsity of aircraft arrivals into a specific flow, it is difficult to obtain meaningful 
and verifiable distributions for aircraft arrival times for individual fiows. Therefore, the aircraft 
arrival process is analyzed over the entire airspace region (i.e. the center level) and proportioned 
amongst the flows. Because the number of aircraft entering the airspace region is a function of the 
time of the day, the objective of this section is to apply the test to determine if the aircraft arrival 
process can be modeled as a non-homogeneous exponential process. Additionally, to strengthen the 
modeling process, similar hypothesis testing is applied to the data to check for a Poisson counting 
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process. While the Poisson counting process implies an exponential inter-arrival distribution, and 
vis versa, testing of arrival counting process yields slightly different results. However, over the vast 
majority of the day, hypothesis testing reveals that both airspace models are appropriate. 

For an arbitrary homogeneous Poisson counting process, N, defined by parameter A 
[events/minute], the probability distribution is given by |26| : 

P(N(i)-N(t-.) = fc) = i^e"^^ (2) 

The associated cumulative distribution for the random variable I representing event inter-arrival 
times, is given by: 

P(I < a) = 1 - e-^^ (3) 

To drive the modeling process, each day is divided into 15 minute increments, during which 
time, the aircraft arrival process is assumed to be static and follow exponential inter-arrival times 
or a Poisson process. The goodness-of-fit test is applied to all 15 minute time intervals for both 
the exponential inter-arrival time distributions, as well as the Poisson counting process. The null 
hypothesis declares that the inter-arrival times into the sector are exponential, with rate ; or in 
the case of the counting process, that the distribution follows a Poisson process with parameter X^t. 
The parameter, X^t, is the average number of aircraft arrivals during the j*'' time interval of length 
T, over all similar days of the week. In this case, r — 15min. For a x^-test result of 1, the null 
hypothesis is rejected, implying that it is unlikely that the inter-arrival times follows an exponential 
distribution (or Poisson counting process) with rate A-' (A^V); otherwise, the hypothesis cannot be 
rejected. The a significance level was set to 0.05, which is a standard value for this test. For the 
exponential distribution, the bin-size for x^-test is selected to be 2 seconds. 



Figure 11(a) illustrates the traffic variation over the day according to the best-fit 15 minute 
parameter A^r for the days Wednesday and Sunday. Over the vast majority of the 24 hour time span 
considered, for Wednesdays (although the same results hold for all other days), the null hypothesis 
is not rejected for either the exponential distribution nor the Poisson process, as illustrated in 
Figure ll(b)| That is, the arrival process for each 15 minute interval follows a best-fit exponential 



distribution. The differences between the values of X^t for Wednesdays and Sundays indicate daily 
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changes in air traffic demand. 

As an example, for the time period between 15:15 and 15:30 Zulu time, over all Wednesday of 
the week, comparison the distribution of aircraft inter-arrival times into the sector, with the best-fit 



exponential distribution, shows a close fit, as illustrated in Fig. 11(c) 

Aircraft arrival process 
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> 

^ 80 
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Time of the day (hours, ZULU time) 

(a) Evolution of the rate parameter (Zulu). 
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(b) t^st results for Wednesdays. 
Distribtuion of Aircraft inter-arrival times ( 15:15 to 15:30 Z) 




■ Data (Wed) 

■ Exponential CDF, 1 11 =14.0027 
1116 inter-arrival samples 

test = 0,p = 0.16731 



40 60 80 

Inter-arrival time Ifk (sec) 



(c) Example emprical and modeled expoential distributions. 
Fig. 11 Aircraft arrival process into the center. 
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The aircraft inter-arrival process into the center is therefore modeled as a non-homogeneous 
poisson process over the complete day. with static rate parameters for any 15 minute interval. The 
rate parameter changes with the time of the day, and day of the week. Additionally, it is likely 
that seasonal changes exist. However for the 4 months considered, any long-term variations over 
same days of the week are small. While, the test does not validate that aircraft arrive with 
an exponential distribution, the test does verify a specific property, enabling modeling of aircraft 
inter-arrival times into the airspace as exponential. Verifying that the air traffic arrival process into 
an en route airspace can be modeled as an exponential distribution supports an assumption found 
extensively throughout the literature (e.g. see |I20I [?fl [25] V 

Similar modeling statements concerning the arrival process into each flow is not possible, as the 
number of aircraft arrivals per flow per 15 minute period is low, and thus fails to yield strong results 
for testing. For modeling purposes, we assume that the flows are independent, and so are their 
arrival processes. The arrival process for the greater airspace is proportioned to according to traffic 
volume along each flow. The rate parameter for the flow Ti during time period j is — tt.^A-', 
where tt^ is given by 

30 

-I = (4) 

d=l 

with ,- being the count of aircraft entering the center during the time period j, on data day d, 
and identifled flying on flow J-";. 

D. Outlier model 

Outliers are aircraft that are not clustered into flows, and they represent approximately 20% 
of the traffic. To model outlier trajectories, density maps of aircraft positions are generated. The 
airspace is discretized into 1 NM squares and 1000 ft high cubes. By spatially interpolating all 
aircraft trajectories, the total number of flights passing through a particular cube is determined. To 
generate a probability distribution, the aircraft counts in each cube are normalized between and 
1 over the entire airspace. Figure [12] presents density plots for the location of outliers at FL 350 
and FL 360. In the flgure, black indicates a high density of aircraft, and white a low density. With 
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exception to the location where major flows pass, the distribution of outliers over the airspace at a 
given flight level is relatively uniform. The spatially probability distribution of outliers deflnes an 
outlier probability model, fo{P)- 

Drawing an analogy with the Kalman fllter, the clusters can be interpreted to be the main 
signal, and the outliers as a noise injected into the signal. 

Outllflrs density fot FL 350 Outliers d en sity for FL seo 




x-axls, NM x-asds, NM 

Fig. 12 Density of outliers at FL 350 and FL 360. 

E. Flow and outlier reconstruction 

Over the last two sections, clustering and outlier analysis provided models for the spatial dis- 
tribution of aircraft through airspace. Furthermore, the vertical and lateral distribution of aircraft 
positions at each window provides greater granularity in the spatial distribution of aircraft. Coupled 
with a model for aircraft speeds, and aircraft inter-arrival times into the airspace, an air traffic model 
is generated. The beneflt of such a model is that it can be parameterized to test 'what-if ' scenarios 
without reconstruction of the traflic model. That is, adjustments to the model properties: flow 
rates, mean aircraft speeds, trajectory windows, or flow trajectories can be made quickly, thereby 
deflning a parameterized generative model. The information provided by parameterized generative 
model is then used to construct the proximity maps proposed in this paper. 
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IV. Generation of aircraft proximity maps 

Modeling airspace via aircraft flows and outliers enables the generation of the presence, conflict, 
and outlier proximity maps, whose computation is now presented. These maps use the information 
generated by the parameterized generative model detailed in the previous sections. We focus first 
on the presence and conflict maps, which consider only aircraft from modeled flows. Then, we will 
introduce outlier proximity maps. 



A. Definition 

Consider the East-North-Up frame TZq = (oo, bo , co ) deflned in Section |II| and fixed with respect 
to the airspace A; and consider any point P — {xp,yp,yp) of the airspace. Given a model of air 
traffic parameterized, the objective is to compute maps based on the probability Pi(P) that at least 
one aircraft flies in the proximity of a point P and the probability P2(P) that at least two aircraft 
from different flows fly in the proximity of a point P. These probabilities will help to generate the 
presence maps (related to Pi(P)) and the conflict map (related to P2(P)). 

The probability Pi(P) can be interpreted as the probability of presence of aircraft, and P2(^') 
can be related to probability of conflict. The shape of the proximity volume is arbitrary but we 
consider flrst the proximity volume as the full cylinder Vc{P) centered on P, and deflned as follows: 

Vc{P) = {{x,v,z)€A\{x-Xpf + {y~ypf<al and {z~Zpf<bl}, (5) 



with ap,6p > 0, as illustrated in Fig. 13 Choosing = 2.5 NM, hp = 1000ft = 0.165NM reflects 
the notion of conflict between aircraft: if two aircraft are simultaneously located inside the cylinder, 
they are in conflict. For any point P, the probability Pi(P) (resp. P2(P)) is then the probability 
that at least one (resp. two) aircraft from different flows are located inside Vc{P), at any given 
time. However, considering the proximity volume VciP) given in Eq. ^ is computationally too 
expensive to generate proximity maps in real-time, since the flows are not only a single track but a 
3D distribution. The cylinder is approximated by its circumscribing rectangular box. The resulting 
values of the probabilities Pi(P) and P2(^) are then upper bounds on the corresponding values if 
Vc{P) had been considered instead. 

For a given airspace geometry and flow conflguration, two three-dimensional proximity maps, 
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Fig. 13 Proximity volume 'Pc(f ) for a given point P of the airspace. 

Ail and A^2, are generated according to the values Pi{P) and P2{P) associated with all points 
P in the airspace. While the map A^i is based on the global presence of aircraft from a flow in 
the airspace, the map M2 highlights the regions of possible conflict between aircraft from different 
flows, which are the most sensitive regions for air traffic controllers. 

B. Inter-aircraft distance and "residual distance" of an aircraft 

From Section [TTj a flow J-i is characterized by, among other parameters, four probability density 
functions : f^k and f^k for the lateral and vertical distributions of aircraft at the window W*^; fATi 
for the inter-arrival time distribution between two consecutive aircraft (for clarity, the index of the 
day j is omitted in this section); fy. for the aircraft speed distribution. To create spatial proximity 
maps, the inter-aircraft distance, i.e. the distance between two consecutive aircraft flying within the 
same flow J-i, must be introduced, instead of inter-arrival time. It will result a complete knowledge 
of spatial distribution of aircraft within the flow J^i, along the vectors (a|, fef, c|). 

Denote {AC^, k E N*} the sequence of aircraft flying in flow Ti. The longitudinal distance Adi 
between two consecutive aircraft of is deflned as the distance between aircraft projection on the 
track (or centroid) 71 (see definition in Section |ll] and depicted in Fig. 14(a) I of two consecutive 



aircraft {AC^ , AC^^'^). The inter-aircraft distance is a random variable such that Adi — ViAti, and 
has associated probability density function fADi ■ Assuming that the two random variables Vi and 
Ati are independent, fADi can be determined from the knowledge of /y^ and /at;- Let Ad™- be the 
mean inter-aircraft distance. 

The notion of "residual distance" of the aircraft is now introduced, as it is a key to the construc- 
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Ad- 



AC! 



(a) Inter-arrival distance 




(b) Residual distance 



Fig. 14 Definitions of the inter-aircraft distance and residual distance. 



tion of the proximity maps. The notion of "residual distance" for aircraft is similar to the concept 
of "residual life" (or "excess life") in renewal theory (see for further details), and it is defined as 
follows: consider any point A located on the track % of fiow J^i (see Fig. |14(b)] ). Observing traffic 
at time to — (we assume that the process starts at time t — —00), let AC^ be the next /c*'' aircraft 
to cross A at any time t > to. The "residual distance" of the aircraft AC^ is the distance between 
the projection of the current position of aircraft AC^ onto % and A. Especially, rj is the distance 
between the current position of aircraft ACf and A. The variable rl is called the "first" residual 
distance because it is related to the first aircraft to cross A. From the PDF of the residual 
distance rf of the aircraft ACf , fjik{r^), can be derived from the PDF of the inter-aircraft distance, 
fADii^di), and is given by the following expression: 



Ad: 



{l-FAoArl)), 



(6) 
(7) 



where Fad. is the cumulative distribution function of the inter-aircraft distance associated to /ad; , 
and * represents the convolution operator. 

To illustrate the notion of "residual distance", consider a straight flow and assume the same 
constant velocity for all aircraft within the flow, given by the maximum velocity estimate of the 
speed distribution. From Section |TT] it is assumed that the inter-arrival time distribution /at; is 
an exponential distribution, characterized by the arrival rate Aj. The PDF /ad; is then also an 
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exponential distribution, characterized by the mean inter-aircraft distance Ad^: 

-^e ^'T"^"' when Ad, > 
fADAAd,) = { ' (8) 
when Ad,; < 0. 

From Eq. ([6|, the probabiHty density function of the "first" residual distance, /^i, can then be 
written as: 

{1 g-si"'-' ^Yien r}>0 
(9) 
when rj < 0, 

which is similar to /ad; for the considered example. For any fc*'* aircraft AC^, the probability 
density function of its residual distance , fj^k, can be then derived from Eqs. Q and ([9|. 

C. Intersection volume 

As mentioned above, considering a cylinder as the proximity volume requires too many compu- 
tations to generate results for a realistic region of the airspace in real-time (i.e. in a few seconds). 
Therefore, the proximity cylinder Vc{P) is approximated by a circumscribing box, whose geometry 
depends of the box considered. Indeed, given any fc*'' box of flow Fi and any point P ^ A. 



the corresponding proximity box V^{P) is defined such that (see Fig. 15 1: 

• its center is P, 

• its dimensions are x x hp (5NMx5NMx0.165NM), 

• it is aligned with the centroid of B^ . 

To compute the probabilities Pi(P) and P2(P), the first step is to determine the intersection 
volume Vf (P), defined to be the intersection of any box with P*^(P), i.e. Vf (P) — B^ H Pf (P) as 



depicted in Fig. 15 Considering both the dimensions of the proximity box (5NMx5NMx0.165NM), 
and the fact that aircraft trajectories are slanted just a few degrees (due to the physical constraints 
of the aircraft fiight dynamics) , it is approximated for computational simplicity that the considered 
box B^ is horizontal around V^{P). 

Consider the frame 7?.^' — (a|,6|,c|) corresponding to the box B^. with P as the geometric 
center. Since TV^ is aligned with the box B^. it is also aligned with Vf (P). The bound of the 
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Fig. 15 Definition of tlie intersection volume Vf (P). 

intersection volume V^{P) is then defined by its entry and exit coordinates along the three axes: A 
and A along a|, _B and B along 6|, and C_ and C along c|. These values can be easily determined 
from the equations of the proximity volume V^{P) and of the box B^. 



D. Presence maps 

Presence maps consider the probability that one or more aircraft is located near a point P, i.e. 
inside the considered proximity volume. Given any point P G consider first a single box , and 
the proximity box Pf (P) . Let Pi(P, B^) be the probability that at least one aircraft ACi from box 
be inside V^{P). Since the three probability distributions for the presence of an aircraft along 
each axis are assumed to be independent (see Section|ll|. the probability Pi(P, Sf) can be written 

Pi(P,Sf) = P(3 AC, e P*^(P)) = P(3 AC, e [AA]).P(3 AC e [BP]).P(3 AC e [CC]). (10) 

So far, only the lateral and vertical PDFs at the windows (i.e. at the entrance of the boxes B^) 
have been defined and determined in Section jllj f^k and f^k . Since the intersection volume (P) 
is not necessarily near , the aircraft spatial distributions within the fiow J-i must be extended 
around the points Therefore, a linear interpolation of the lateral (resp. vertical) distributions 

f^k and f^k+i (resp. f^k and f^k+i) between the windows and W*^^^ is used to get fa^/A (resp. 
fck/c)^ ^-^^ lateral (resp. vertical) aircraft distribution at the point A (resp. C). The linear 
interpolation process approximates the "smooth" variation of aircraft presence density within the 
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flow, preserving the nature of a probability density function. 



Any term in the right-hand side of Eq. ( 10 1 can be computed as a function of the probabiUty 



density functions previously defined, /^fc/yi, /c?f/ci Ir^ '■ 



Pi3 AC, e[AA])= / f,.^A^)dx, 

J A 

P{3 AC, e[CC])^ r f,.,c{z)dz, 
Jc 



P(3 AC, G [BB]) 



BB 



Ir} (y) dy, 



(11) 

(12) 
(13) 



where BB = 



BB 



The probability Pi(P,Sf) can then be determined from Eqs. (|lO|-fL3l. 
The previous definitions and equations can naturally be extended to the li boxes of the fiow 
J^i. It leads to introduce the probability, Pi(P, J^i), that there is at least one aircraft ACi from fiow 



J^i flying at proximity of the point P. The probability is computed from Eqs. ( 10 1-( 13 1 taking the 
li boxes into account. 

Consider now the entire airspace A modeled with N independent flows. The presence probability 
of at least one aircraft from a flow about any point P, Pi(P), can then be written 



N 



N 



Pi(p) = 1 - p(yi e [1,7V], $ AC, e r,{p)) ^i^Y[Pi$ AC, e r,iP)) - 1 - n(i - Pi(^,-^^))- 

(14) 



i=l 



E. Conflict maps 

Conflict maps determine the probability that at least two aircraft from different flows are simul- 
taneously in the same region of the airspace. More formally, the probability P2(P), that at least 
two aircraft from different ffows are in the proximity volume surrounding P, is written : 

P2(P) = P(3(z,j) e [l,N]^,t^j, s.t. {AC,eV,{P)) k {AC,eVj{P)). 



Using the probability Pi(P) deflned in Section IV D P2(P) can be written 



P2(P) = Pi(P) - P(3i G [1,N], (AC, e V,iP)) k [ACj i P,(P),Vi ^ i)) 

N N 

=Pi(p)-^Pi(p,j-,) n (i-Pi(^,-^,))- 
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F. Outlier proximity maps 

Presence and conflict maps consider only aircraft that belong to flows. For a complete modeling 
of the traffic, outliers - non clustered aircraft - should be considered because they account for 
20% of the overall traffic and may have an important influence on the dominant traffic. The 
spatial distribution of the outliers, /o(-P), is given by the parameterized generative model presented 
in Section This leads to the introduction of outlier proximity maps, denoted A4o- Outliers 
proximity maps estimate the influence of outliers on the modeled flows that represent the majority 
of traffic inside the airspace under consideration. Therefore the probability Po{P), that one aircraft 
from the flows and one outlier are simultaneously in the same neighborhood of any point P of the 
airspace, is determined. Only one outlier is considered, since the probability of two outliers to be in 
proximity too each other is negligible. Using a methodology similar to the one used in Sections |IVD| 



and IV E Po(P) is approximated by 



1 

ipbp J J Jvo[P) 



Po(P)=Pi(P).^ /// fo{P')dT, 



where Vo{P) is a parallelogram whose dimensions are ap x ap x 2bp, and whose center is P. The 
spatial distribution of the outliers, fo{P), can be designed (i.e. simulated) to evaluate the influence 
of the presence of rogue aircraft (considered as external disturbances) on a "nominal" air traffic 
consisting of the flows in a given region of the airspace: this approach extends the works in |12l 113] 
from individual aircraft to entire aircraft flows. 

V. Generation of proximity maps for Cleveland center ZOB 

This section combines both historical flow parameters determined from the clustering algorithm 
presented in Sections [IT] and |III| and the methods to construct proximity maps given in Section |IV[ 
The maps use the flows characteristics computed for May and the time period between 2:00p.m. 
and 2:15p.m. local time. During this time period, the expected number of entering aircraft is 65. 
For easier visualization of the results, proximity maps have been computed for the horizontal slice 
of airspace between FL 310 and FL 350. In this region, the mean inter-aircraft distances for the 10 
busiest flows are given in Table [T| and the centroids of the flows are plotted in red on Figs. 16(a)| 



and [T7(a)] Fi gures [l6| and [T7| illustrate the presence and conflicts maps, respectively. A 3D view 
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and a 2D view are presented for each flight level. The outliers proximity maps for the considered 
center are given in Fig. [18] which highlights the regions of the airspace that may be sensitive to the 
presence of an outlier. These maps provide an intuition of the importance of each flow, how they 
interact with each other, and how much influence an outlier can have on the main traffic. Since 
these maps provide a depiction of the probability of proximity of aircraft, they may be used in 
combination with a forecasted weather map, in order to estimate the regions where bad weather 
may have strong impact on the predicted traffic (as illustrated in [21). The maps can also be used 
as a tool for traffic flow managers to redesign routes or adjust flow rates in order to allow increased 
throughput while keeping a low level of proximity probability. 



Flow number, i 
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Table 1 Mean inter-aircraft distance, Cleveland center ZOB between 2:00pm-2:15pm 

VI. Conclusion 

This paper is concerned with computing airspace statistical complexity metrics using a parame- 
terized generative model driven by historical flight data. Trajectories are clustered into major flows 
that account for 80% of the traffic and the remaining 20% are identified as outliers. These flows and 
outliers constitute a structure of air traffic, used as a parameterized generative model. Within the 
flows, aircraft spatial and temporal probability distributions are computed. Analysis of the airspace 
demonstrates that the aircraft arrival process into the center can be modeled as time- varying Poisson 
process. Using the structure determined by the flow and outliers model, three-dimensional maps 
are constructed to predict and visualize the probability of presence and probability of proximity of 
aircraft in the airspace. 

The proximity maps are information displays that are relevant to traffic flow management. Each 
map (presence, conflict, or outlier proximity map) identifies and highlights "critical" regions of the 
airspace that require signiflcant monitoring. Further, because the maps are based on parameterized 
historical traffic models, it is possible to consider "what-if" scenarios, in which parameters (flow 
rates and distributions) are adjusted in real-time to consider daily traffic variations. Ultimately, 
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Presence Maps for Celveland Center ZOB 




Presence Map for FL 31 



Distance y, NM 



Distance x, NM 



(a) 3D view 
Presence Map for FL 330 
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(b) 2D view at FL 310 
Presence Map for FL 350 
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(c) 2D view at FL 350 (d) 2D view at FL 380 

Fig. 16 Presence maps of Cleveland center ZOB between 2:00pm-2:15pm. 



this tool provides an information display to aid in the strategic and operational decision making 
required in traffic flow management. 
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Fig. 17 Conflict maps of Cleveland center ZOB between 2:00pm-2:15pm. 
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(a) 2D view at FL 310 (b) 2D view at FL 350 

Fig. 18 Outlier proximity maps of Cleveland center ZOB between 2:00pm-2:15pm. 
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